Solvable dynamics in a system of interacting random tops 
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In this letter a new solvable model of synchronization dynamics is introduced. It consists of a 
system of long range interacting tops with random precession frequencies. The model allows for 
an explicit study of orientational effects in synchronized phenomena. A stability analysis of the 
incoherent solution is performed for different types of orientational disorder. A system with only 
orientational disorder always synchronizes in the absence of external noise. 

One of the main goals in modern statistical physics concerns the understanding of dynamical behavior in complex 
systems formed by a large population of interacting entities. Much theoretical effort has been devoted to the study of 
synchronization dynamics in simple systems composed of interacting units. In those cases there is competition between 
oscillations arising from the natural randomness in the members of a population and macroscopic synchronization of 
the population as a whole. It is widely believed that these models provide a plausible explanation for the existence 
of synchronization phenomena in a large variety of physical systems ranging from biology to economyEI. On the 
other hand, the study of this type of models provides an interesting framework to study relaxational processes and 
off-equilibrium behavior in physical systems driven by external forces. 

A simple model which describes the emergence of synchronization phenomena in a population of phase oscillators 
was proposed many years ago by KuramotoB. In this model each oscillator is characterized by a random natural 
frequency and interacts ferromagnetically with the rest by a mean-field interaction. The ferromagnetic ordering of the 
system competes with the independent random oscillation of each oscillator. Despite extensive studies in the past, 
there are still some open issues such as aijull characterization of the stationary synchronized states, tharstudy of the 
dynamics in the absence of external noisea as well as the onset of synchronization in the critical regionQi3. 

The aim of this letter is to introduce a new model which yields insight on the mechanisms involved in synchronization 
phenomena. The fundamental new feature of this model is that it explicitly introduces the role of orientational 
degrees of freedom in the synchronization dynamics. This feature is not cajstured by the Kuramoto model and can be 
interesting from the viewpoint of biological systems and also neurosciencea where orientational effects are certainly 
important. 

The model consists of a system of N tops, each one characterized by a random natural precession vector cJ^, 
interacting ferromagnetically in a mean-field way. The tops are specified by a three components unit vector Xi 
(i = 1, N) and obey the following dynamics, 

dxi dH ^ ^ ^ 

where the Coi are random quenched precession vectors, TL is the Hamiltonian of the system and fji is an external white 
noise with zero mean and correlation {'f]i(t)ffj{t')) = 6TSijS{t — t'). The Hamiltonian TL is given by 

rL = -^^XiXj ~^^l^xj (2) 

i<'j i 

where the fii are Lagrange multipliers introduced in order to ensure the local constraint \xi\ — 1 at all times. 
Note that in the dynamics there is competition between two opposite effects: a ferromagnetic interaction which 
tries to align the tops in the same direction and a natural precession of the tops around random quenched di- 
rections which drives the system to an incoherent state. To analyze the previous dynamics it is convenient to 
introduce polar coordinates for the tops Xi = (sin 0^ cos 0.;, sin^i sin (/>i, cos 0i) as well as for the random precessions 
LUi = uJi{sm fii COS Xi, sin Hi sin Xi, COS fii). It is not difficult to check that dynamical equation (|^) can be written in the 
following way. 



dt 



= —K{C sin 6*2 -f A cos 9i cos (pi + B cos 6*2 sin cpi) + uji sin jii sin(Ai — (pi) + ai (3) 



sin^i = —K{Asin(l)i ~ B cos (pi) + w,;(— sin^^ cos Oi cos(Ai — (pi) + cos fii sin 9i) + (3i (4) 
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where NA = J^j cos sin 6j , NB = J^j sin (f'j sin Oj , NC = J^j <^os 9j , and the ai , /3i are Gaussian noises given by, 

(a,{t)) =Tcot9f, {a,{t)aj{t')) = 2TS,,5{t - 1') (5) 

(A(t))=0; mt)l3j{t'))^2T5,,5{t-t') . (6) 

Note that the mean vahie of the noise a is not zero but depends on the polar angle 8. This is expected because 
the angle 6 lies always within an interval (0, tt) of length smaller than 2tt and = 0, tt are returning points of the 
dynamics. Solving the previous dynamical equations seems at first glance a rather difficult task. Here we will follow 
a powerful approach recently introduced for the study of the Kuramoto modelEI by considering an appropriate set of 
moments which is invariant under the symmetry of the original dynamical equations (|l|). Before studying the most 
general case, and for sake of simplicity, we will first consider the case where there are no frequencies, i.e. Ui = 0. 

The non-disordered case toi = 0. It is easy to observe from equation that dynamics of the model is invariant under 
the group of spatial rotations whose generators are and and the eigenfunctions are the spherical harmonics. 
The most natural set of moments which can close the dynamics is, 

1 ^ 

MUt) = j;^Y.YiMt),Ut)) (7) 

i=l 

where the Yira{0, (j)) = CimPi^{d) exp{im(j)) are the spherical harmonics with Cim — ((2? + — TO)!/47r(Z + m)l)i as 
normalization constants. Pj"* are the associated Legendre polynomials where —l<m<l. 

It is not difficult to write the equation of motion for the Mi^lt). Using the definition of spherical harmonics, simple 
recursion formulas for the Legendre polynomials as well as the Gaussian representation of the noise, it is possible to 
show that the moments obey the following closure equation, 



at 



,m+l ^l.m 

+ iA + iS)(e,,™Af/_i,„,_i + fLmMl+i,m-l)) - Tl{l + l)Ml.,n (8) 



where the a, 6, c, d, e, / are Clebsch-Gordan-like coefficients given by. 



_ + 1) j {l-m){l~m-l) _l_ l {l + m+l){l + m + 2) 

2 V 4F-1 ' '^'•"'^2V (2/ + l)(2/ + 3) ^ ' 



and ei^m — ci_-m, fi,m = di^—m- The time dependent parameters A, B, C in (^) are given by A+iB = V8^Mi,i , C = 
a/ 47r/3Mi_o- Note that the dynamical equation (||) is invariant under the transformation m — m if the moments 
satisfy the relation Mi^-m = (—1)™ ^^*m which is indeed the relation satisfied by the spherical harmonics F;,™. The 
recursion relations eq.(^) show explicitly that the dynamics has been closed. TJbe special case K = corresponds 
to the random walk on a spherical surface, a case well known in the literaturaa. Then eq.(^) trivially reduces to 
Mi^rn — ~-Tl(l + l)Mi^rn which shows that all moments decay exponentially fast to zero (except Mq.o which is a 
constant of motion equal to I/V^tt). Eq.(H) can be expressed in a more intuitive form introducing an appropriate 
generating function. To see this, let us define 

00 i 1 "'^ 

gt{e, 4>)^Y.Y. MunmiMO, = ]^ E '^(co^^^ - COS 0)5(0, - </.) (11) 

1=0 m=-l j = l 

which is nothing else than the probability density to find a top with a given solid angle il = (0, 4>) on the unit sphere. 
The last equality in eq. (pi]) comes from the closure condition of the spherical harmonics. Using eq.(|^) it is not difficult 
to derive the following equation for gt: ^ — ~div{gtv{6, <j))) + T'V^gt where v — {ve, v^) is a two component velocity 
field given by, vg — iirr(sinOcos6'cos($ — </)) — sin 6* cos 0); — ifr sin O sin(<I> — 0). In polar coordinates the 
equation for gt is: 
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The time dependent parameters r, 0, $ which appear in the velocity field are self-consistently computed from the 
probability distribution gt. They are given hy A = r sin(0) cos($); i? = r sin(0) sin($); C = cos(O) where A^B^C 
have been introduced before, r = ^ A^ + B"^ + is the synchronization parameter and measures the degree of 
coherence of the tops. Equation (^2|) seems hardly manageable but still some results can be inferred, in particular the 
nature of the stationary solutions. It is easy to check that these are solutions of the Boltzmann type, i.e. equilibrium 
solutions of the Hamiltonian eq.(^. This is expected since the model, in the absence of random precessions, is purely 
relaxational. We will not extend further on the connections between the moment approach shown here and the 
probability density formalism. The equation (^2|) is formally identical to that derived for the Kuramoto model with 
the only difference that its velocity field has now two components. Certainly it can be generalized for 0(n) vector 
mean- field models with any number of components (n=2 is the Kuramoto model and n=3 the present case). Since 
our main interest in this note is to analyze the model with random precession frequencies, we will come back to the 
moment approach which is more appropriate and suited to investigate the dynamical phase diagram of the disordered 
model. 

The disordered case uji ^ 0. In the presence of random precessions the model is not purely relaxational since there 
are external driving random forces. In this case we expect the emergence of a rich dynamical behavior due to the 
competition between the ordering ferromagnetic interaction and the random natural precessions of the tops. For the 
sake of simplicity and in order to investigate the effect of orientational disorder, we will consider here the case in which 
the precession angular velocities have the same magnitude (i.e. LOi = lo), but point in different directions in space. 
In this case the disorder is specified by a probability distribution p{fi,X). It is easy to generalize the definition of 
the mornents (|^) to include the presence of quenched disorder. Now the moments are characterized by four quantum 
numberso, two of them appearing as a consequence of the disorder, 

1 ^ 

= nT. A.) . (13) 

After some algebra, the closure equations for the new set of moments read 



K{C [ai^MfXm - h^Mf^lJ -{A- ^B){curaMfX„^^, + d,™Mf4,„+i) + 
{A + zi3)(e,,„Mn,,„_, + /,™Mf/,,,„_ - Tl(l + mfl + Y^'^ZmC''''' (14) 

where the expressions for a, 6, c, d, e, / have been given in eq.(^ 10) and the other coefficients are 



.P,q _ l {l-m + l){l + m){p-q)ip-q-l) _ j {p + q + l){p - q + 1) 

'■™"V V_i ' (2p+l)(2p + 3) ^^^^ 



_ + m){p + q + 2){p + g + 1) . _ j [p - q){p + q) 

'■™"V (2p+i)(2p+3) ' ^^-^'-y V_i ^^^^ 

and Cf^ — A^'Zmi ^Tm ~ ^f'-m- "^^^ parameters A,B,C have been already defined after eqs.(||,^). In terms 
of the new set of moments they are given by, A + iB ~ 47r-y/2/3M°'^; C — 47r-\/l/3Af^'g . Note that the set of 
moments Mqq — -^j^ J dilYp^q(n)p{n) (where 51 — (/i, A) is the solid angle and the integration is over all the 
unit sphere) are constants of the motion. Let us mention that also in this case a probability distribution, i.e. a 
generating function for all the moments, can be defined like in the non-disordered case. We will not extend on these 
considerations and instead we will focus on the novel properties of the model. For this purpose let us consider the case 
of disorder distributions with axial symmetry around the z-axis (the following considerations can be easily extended 
to the more general casell3). This means that p{n) = p{l^)- In this case it is possible to study the region in the 
plane oj — uj/T,K — K/T where the incoherent solution r = is unstable. Our calculation follows essentially the 
equivalent one for the Kuramoto modelliilllj. To this end we expand the moments around the incoherent solution 
Mf"^ = -^=5ifi5mfiM^'^Q + er]f'^e'}q){at) where e is a smaU parameter. Substituting this result in equation ( |l^ ) 
and performing a linear stability analysis we find a set of two linear equations (uncoupled from the rest of modes) 
involving the two fundamental modes 771^0(^^)1 '7i,i(fl)- These modes are defined by, r]f'^ = J dflr]i^mi^)Yp^q(Q) . After 
some calculations (details will be shown elsewheret2l) the condition for the linear stability of the incoherent solution 
is determined by the roots of the equation Ai(7)A2(7) = where. 
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\ih)^l^-—l^+Co^^-—Cj^cos^ifi) (17) 

2K K 

A2(7) =7^- — 72+(I,2^-— w2sin2(^) , ^ = a + 2 (18) 



and A{fi) — 2tt sm{fi)p{fi)A{fi). These are two cubic equations and each one yields three roots (one of them is real 
and the other two are complex conjugates). When the real part of one of these six roots becomes larger than 2 the 
incoherent solution is unstable. The equation in (|l^Jl^) that determines the stability is the one where the independent 
7 term is the largest (see the examples below). The boundaries of the region where the incoherent solution is 
unstable indicate the dynamical transition lines. While the stability boundaries will depend on the particular disorder 
distribution, some general results can still be inferred. In particular, in the absence of noise T = 0, there is no critical 
value of K below which the incoherent solution becomes linearly stable. This result is also foupd in the Kuramoto 
model with disorder pdjstributions with vanishing probability of oscillators with zero frequencytiJ, for instance, the 
bimodal distributionEj. This last case has been shown to display a very rich dynamical behavior also shared by 
the present model. Here we are going to consider three cases of disorder distributions which cover a large variety of 
physical situations: a) precession vectors lying in the z-axis and randomly pointing in opposite directions, i.e. p{fi) = 
except for /j, = 0,7r; b) precession vectors uniformly distributed and lying in the x-y plane, i.e. p{fj,) — c) 
precession vectors isotropically distributed on the sphere, i.e. p(/i) = 

Case a) is particularly interesting despite its simplicity. Half of the tops precess in one sense, the other half do in 
the other sense, around the z-axis. In this case cos^di) = 1 and the stability condition for the incoherent solution 
is determined by Ai(7) which yields < 3. For K > 3 the incoherent solution is always unstable whatever the 
value of w. A quick inspection on the definition of the moments cq.(^3|) reveals that they reduce to two families 
Hi^m = ^, Gi^m = The dynamical equations for the moments H and G are very similar to those obtained 

for the Kuramoto model in case of a bimodal distribution^, the crucial difference being that the two families of 
moments H and G are coupled each other with a term proportional to the quantum number m. A trivial stationary 
solution of dynamical equations is then given by Hi^m — ^/.o<^m,Oj Gi^m = 0. This solution coincides with that derived 
previously in the non-disordered case uJi — 0. Interestingly enough we find that the stationary solutions of the model 
(a) are the equilibrium solutions ofptiie Hamiltonian eq. (y) even though there is no effective Hamiltonian which yields 
the original equation of motion (|l])t2l. 

Cases b) and c) display a more rich dynamical behavior. Case b) corresponds to all precession directions lying on 

the x-y plane. In this case sin^(/i) = 1 and the stability criteria is fixed by the roots of A2(7). The incoherent solution 
is linearly stable for if < 3, linearly unstable for K > 6 and for intermediate values of K depending on the value of 
a). More concretely, the incoherent solution is stable for lj^ > {8K — 24)/(6 — K). Case c) is the limit situation (i.e. 
the most orientationally disordered case) with the largest region in the phase diagram where the incoherent solution 
is linearly stable. In this case, cos2(/i) = 1/3 and Ai(7) = A2(7). The incoherent solution is linearly stable for K < 3, 
linearly unstable for K > 9 and stable in the intermediate region 3 < K < 9 ii Cj^ > {12K — 36)/(9 — K). 

The main difference between the present model and the Kuramoto model with bimodal distribution of frequencies 
lies in their phase diagrams. The neutral stability line for the incoherent solution in the phase diagram (K, Co) has a 
vertical parabolic branch at Kc- uj ^ {K~Kc)~^ {Kc =pj3, 9 for cases b) and c) respectively). In the Kuramoto model, 
the neutral stability line is the semirect Kc = 4, a) > lE3. In the region where the incoherent solution is unstable we 
expect the existence of coherent fixed point solutions as well as oscillatory time dependent solutions. To show this we 
have numerically solved eq.(|lj) for a finite set of moments using a simple second order Euler algorithm. A maximum 
value of both angular momenta I an p equal to 10 is enough for the calculation (which implies an approximate number 
of 10000 moments). In figure 1 we show the time evolution of the synchronization parameter for model c) in the 
region where the incoherent solution is unstable starting from two different initial conditions. 
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FIG. 1. Synchronization parameter r as a function of time for case c) at K — ui = 50/3 starting from two initial conditions: 
all tops pointing in the same directions (continuous line) and the uncoherent state (dotted line). 

To summarize, we have introduced an solvable model of interacting random tops. The model explicitly introduces 
orientational effects in the synchronization dynamics, a feature which is not present in the Kuramoto model. The 
model has been analytically studied in the orientationally disordered case introducing a family of moments (depending 
on four quantum numbers, two of them arising from the randomness in the system) which can be exactly closed due 
to the symmetry invariances of the original dynamical equations (|l|). The moment formalism is simple and suited 
for theoretical analysis and this has been explicitly checked in the study of the stability properties of the incoherent 
solution. We find that the system with only orientational disorder always synchronizes in the absence of external 
noise. The present model (n=3) altogether with the Kuramoto model (n=2) provide simple examples of analytically 
tractable models with non relaxational dynamics. There still remain some open issues like the study of the effect of 
frequency dispersion (i.e. the Ui different) as well as the noise free (T = 0) dynamics in this model where the higher 
modes of the incoherent solution are neutrally stable. 
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